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Abstract 

> 

We address the problem of simulating pair-interaction Hamiltonians in n node quan- 
tum networks where the subsystems have arbitrary, possibly different, dimensions. We 
show that any pair-interaction can be used to simulate any other by applying sequences 
of appropriate local control sequences. Efficient schemes for decoupling and time rever- 
sal can be constructed from orthogonal arrays. Conditions on time optimal simulation 
are formulated in terms of spectral majorization of matrices characterizing the cou- 
pling parameters. Moreover, we consider a specific system of n harmonic oscillators 
with bilinear interaction. In this case, decoupling can efficiently be achieved using the 
combinatorial concept of difference schemes. For this type of interactions we present 
optimal schemes for inversion. 
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1 Introduction 

The conjecture that quantum computers might be able to simulate the time evolution 
of quantum systems better than classical computers has already been stated in ||. 
Various schemes for constructing gate sequences which simulate the unitary evolution 
corresponding to a given Hamiltonian have been suggested (see e.g. [18, 23|). More 



recently, a quite different approach to this problem has become popular: considering 
the quantum computer as a quantum system with Hamiltonian evolution as well, a 
simulation is a sequence of control operations acting on the quantum computer in such 
a way that the net effect is a time evolution analogous to the evolution of the system 
we want to simulate || [12], [j], |2^]. In the setting described in the following the 
simulation problem can be stated in a control-theoretical fashion. Assume that the 
total quantum system is a quantum network, i. e., a system consisting of n subsystems. 
The Hilbert space is the tensor product 

H:=Hi®H 2 ®---®H n , 
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where[] Hj := C d for j = 1, . . . , n and d € N. Let {<r a | a = 1, . . . , d 2 — 1 =: m} be 
an orthogonal basis of the R-vector space su(d) of traceless Hermitian operators on 
C d . We assume that both the Hamiltonian to be simulated and the Hamiltonian H 
of the system are sum of pair-interactions between nodes and free evolutions on each 
individual node. Hence the system Hamiltonian is given by 

kl;af3 k;a 

where J is a real symmetric mn x mn-matrix and r is an mn-dimensional real vector. 
Furthermore, we assume that the only transformations which can be implemented 
directly by external control interactions are local transformations of the form 

U := Vi ® V 2 <g> • • • <g> V n , (2) 

where each Vj for j = 1, . . . , n is an element of the special unitary group SU(d). As- 
suming that the implementations of Vj are fast compared to the natural evolution given 
by H ("fast control limit"), the time evolution according to U^HU can be simulated 
by alternating the natural time evolution with implementations of and U. Here we 
make use of the identity 

exp(— iHt)U = exp(-iU^HUt) . 

Concatenating the unitary transformations exp(— if/J HUjTj) we obtain an approxima- 
tion of the time evolution corresponding to the "average Hamiltonian" 



if the times r,- are assumed to be sufficiently small (for details see Jq])q. In this sense, 
the set of Hamiltonians which can be simulated with no time overhead is exactly the 
convex span of the set S := {U^HU}, where U is as in eq. (pf) . For a formal definition of 
the notion of time overhead of a simulation see [p7| , |l| . Intuitively, if the only possibility 
to write a Hamiltonian as a positive linear combination of elements of the form U HU^ 
is to do it in such a way that the sum of the coefficients is greater than one, then the 
sum of these coefficients is precisely the time overhead of the simulation. A suitably 
rescaled Hamiltonian is then in the convex span of S. As noted in p7| , the minimal 
overhead for simulating a Hamiltonian H by the physical Hamiltonian H is the smallest 
positive r such that H/t is in the convex span of S. 

Note that convex problems of this kind are closely related to the method for obtain- 
ing pseudopure states by averaging over random unitary transformations [16|. Pseu- 



dopure states are states that can be written as convex combination of the maximally 
mixed state (with density matrix 1/d) and a pure state {i/j)^. Writing a general state 
as p = 1/d + A, where A is the traceless part, we have that p can be transformed into 
the pseudopure state (1 — A) 1/d + X\ip) by averaging over unitary transformations 



1 For reasons of convenience of notation we assume that the subsystems have equal dimensions. However, 
note that our results on universal simulations generalize straightforwardly to arbitrary dimensions. 

2 Note that an approach of this kind is generally accepted for describing Nuclear Magnetic Resonance 
experiments. 
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if and only if A can be transformed (without time overhead) into the traceless opera- 
tor —Xl/d + A\ip){ip\. Determining the optimal signal-to-noise ratio of the attainable 
pseudopure state is hence directly related to determining minimal time overhead of 
simulation schemes. 

For an arbitrary interaction between two qubits, minimization of the time overhead 
has been carried out in [jy]. In case of n qubits and pair-interaction Hamiltonians H 
and H the situation is more complicated: the problem of optimal simulation of H by H 
cannot be reduced to the two-qubit case. The difficulty arising here is that the required 
operations on qubit k for simulating an interaction between qubit k and I might not 
coincide with the required operations on k for simulating another interaction between 
k and m. The fact that the control operations for simulating interactions between 
non-disjoint qubit pairs must be consistent on the common qubit seems to be a highly 
non-trivial combinatorial problem. Some upper and lower bounds for the time overhead 



are given in Jl^, 27]. 



The issue of optimal simulation has consequences for the problem to parallelize 
operations of discrete quantum algorithms. The maximal degree of parallelization is 
constrained since a given interaction may allow to implement concatenations of some 
two-qubit gates even if they act on non-disjoint qubit pairs. Assume for instance that 
a Hamiltonian on three qubits is given by 

H := a z (g> a z (g> 1 + o z ® 1 ® o~ z + 1 (g) o z (g> a z , 

where a z denotes the Pauli matrix. Then, for generic times t > 0, the unitary 
exp(-iHt) is a concatenation of three two-qubit gates, which can be implemented 
"simultaneously" by waiting the time t. But there is no obvious way for such a fast im- 
plementation of exp(-iHt) for negative t provided that H is the system Hamiltonian. 
One can even show that there exist negative values of t such that the implementation 
of exp(-iHt) cannot be performed in principle within a time period of length \t\, since 



H is not able to simulate —H without time overhead [12|. For n qubits with a z (g) o~ z - 
interactions the minimal overhead is known to be at least n — 1. This shows that the 
question of optimal parallelism of sequences of gates is rather sensitive to the form of 
the underlying Hamiltonian. Hence the question of optimal implementation of com- 
plex transformations on an n-partite system cannot be answered alone by resolving the 
network into elementary gates: the question of maximally possible parallelism appears 
already on the control-theoretic level. 

The paper is organized as follows. In Section [2| we construct decoupling schemes 
from orthogonal arrays. Selective decoupling can be achieved by straightforward gener- 
alizations. Therefore each pair-interaction Hamiltonian on an n-partite systems can be 
converted into any other provided that a sufficient set of local unitary control operations 
is available (the so-called transformer groups introduced in [p8[ ). Such a simulation has 
time overhead of 0(n 2 ). Our schemes works even if the dimensions of the n subsystems 
do not agree, thereby solving a problem stated in [poj ], 

A lower bound on the time complexity of mutual simulation of different Hamilto- 
nians is shown which uses the spectra of matrices describing the coupling parameters 
of the interaction. Efficient schemes for switching off all the interactions ( "decoupling 
schemes") are constructed. Using the combinatorial concept of orthogonal arrays de- 
coupling can be achieved if local transformations of a unitary error basis are available. 
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The length of the required sequence of local transformations grows linearly with the 
number of nodes. 

A specific form of interactions is dealt with in Section ||. There we consider harmonic 
oscillators coupled by bilinear terms of creation and annihilation operators. In this case 
rather simple decoupling schemes can be constructed on the basis of difference schemes, 
another concept from combinatorics. Furthermore we discuss the generalization of this 
problem where a given bilinear interaction is used to simulate another interaction of 
this form with different coupling parameters. Finding the time optimal simulation 
is equivalent to the mathematical problem of constructing vectors of minimal length 
with complex entries of modulus one such that their inner products yield certain values. 
This fact is used for deriving a lower bound on the simulation time. In contrast to the 
general system Hamiltonian, the upper bound is O(n). An optimal scheme for time 
inversion is constructed. 

In Section ||| we discuss the connection between the decoupling schemes of Section |2| 
and the schemes presented in [|T7|]. 



2 Simulating Hamiltonians in Networks 
2.1 Selective Decoupling 

A useful tool for simulating Hamiltonians in multi node quantum networks is given 
by decoupling schemes, i.e., sequences of local operations which switch off unwanted 
interactions. First we describe schemes for decoupling all the interactions. Schemes for 
decoupling two interacting quantum systems — mainly in the context of a quantum reg- 
ister which is coupled to a bath — have been derived in |2^]). Decoupling of n-partite 
systems interacting with each other by two-body Hamiltonians has been considered in 
1 17, 24] for the case of qubits. In it was noted that decoupling of bipartite systems 
can be concatenated in such a way that the system is separated into suitable clusters 
of subsystems without coupling between different clusters. The number of operations 
involved in this scheme is 0(n 21og2d ). In this paper we present a construction based 
on orthogonal arrays which uses O(n) operations provided that the dimension d of the 
nodes is a prime power. 

Straightforward generalizations of decoupling schemes allow to switch off all in- 
teractions except the Hamiltonian of a pair of nodes ("selective recoupling") or the 
Hamiltonian of a single node. We have shown in (2^] that any bipartite Hamiltonian 
can simulate any other provided that it consists of non-trivial local Hamiltonians on 
both nodes and a non-trivial coupling. If this criterion is met by all pairs of an n- 
partite Hamiltonian then universal simulation of all pair-interaction Hamiltonians is 
possible^. This result is also true if not every local operation is available, it is sufficient 
that all the elements of a so-called transformer group (a concept introduced in |2q ]) 
can be implemented. The simulation time overhead is of order 0(n 2 ) since there are 
n(n — l)/2 pairs of nodes. 



Rephrasing well-known results about decoupling (e. g. |26|, |29(]) in our language, we 



3 Note that the condition that all the couplings have to be non-trivial is only necessary in the average 
Hamiltonian approach. If higher order terms in the time interval are considered interactions between nodes 
k and I resp. k and m can be used for simulating a coupling between I and m as noted in J^, |20fl . 
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briefly describe decoupling for bipartite systems. Remember that the Hamiltonian on a 
d-dimensional quantum system can be cancelled if an appropriate sequence of unitary 
control operations is applied p8[|. 



Definition 1 (Annihilator) An annihilator A := (17%, n, U2, t%, . . . ,Un,tn) of di- 
mension d and length N is given by unitaries XJ{ 6 SU (d) and relative times T{ > 0, 
YliLi r i = 1 such that 



N 



YjTiUlaUi = 



i=l 



for all a G su(d). An annihilator is called minimal if there is no shorter annihilator. 



A minimal annihilator of dimension d has length d 2 (as already noted in [26] and 
proved in ^8|) and all 7$ are equal. Moreover, the unitaries Ui must form a unitary 
error basis of operators in C dxd , i. e. a collection of d 2 unitaries U that are orthogonal 
with respect to the trace inner product (A\B) := tv(A^ B)/d. Such unitaries can be 
explicitly constructed using nice error bases JIB] , |i~4j ] . 

Being a special case of equation (|l|) we write the general Hamiltonian of a bipartite 
system as follows: 

H = J Q/3 cr Q ® 073 + r 7 1 (g> cr 7 + s s cr 5 (8) 1 . (3) 

a/3 7 <5 

Let £1 = {C/j} and £2 = {^i} be unitary error bases of the respective systems and 
let A denote the set {1, . . . ,d 2 }. By applying the annihilators defined by E\ and £2 
independently on the nodes we can switch off the Hamiltonian, i. e. 

T^i E (u!®Vj)H(U i ®V 3 ) = 0. (4) 

We describe a decoupling scheme on n nodes by n unitary error basis £1, ■ . . ,£ n 
and an n x iV-matrix M = (772.17)1=1,... nj=i... N- This matrix contains elements of A 
specifying the conjugation by the mfj unitary of £{ on a specific node % for a certain 
time interval j. The iV time intervals correspond to the columns and the different 
nodes correspond to different rows. For instance the decoupling scheme corresponding 
to (|j) is given by the array 



11... 1 
1 2 ... N 



2 2 ... 2 
12... N 



and unitary error bases £\ = {U} and £2 = {Vi} corresponding to the first and second 
node, respectively. 

The simplest approach for decoupling is to choose the columns of M as all tuples of 
A n . However, this scheme is not efficient in terms of the number of time intervals and 
pulses since both scale exponentially as d 2n because the sequence has to be repeated 
d 2 times for each added node. 

More efficient schemes can be constructed using the combinatorial structure of 
orthogonal arrays. See ||, ||, |9| for the general theory of orthogonal arrays. Orthogonal 
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arrays have numerous applications e.g., in the design of experiments. Also there are 
connections between orthogonal arrays and mutually orthogonal Latin squares and 
transversal designs (cf. ||, Section VIII]). The following definition takes account of 
the fact that for purposes of decoupling we need a special type of orthogonal arrays, 
namely those of strength t = 2 (cf. |2|, ||, || for the general case). Also the notation 
used is adapted to this situation. 

Definition 2 Let A be a finite alphabet and let n,N £ N. An n x N array M with 
entries from A is an orthogonal array with \ A\ levels, strength t = 2, and index A if and 
only if each pair of elements of A occurs A times in the list ((mkj,mij) \ j = 1, . . . N) for 
1 < k < I < n. We use the notation OA\(n, N) to denote a corresponding orthogonal 
array^ 

The following theorem shows that decoupling in networks of arbitrary dimensions 
can be achieved using pulse sequences obtained by orthogonal arrays. 

Theorem 1 (Decoupling) Let A be the finite alphabet {1,... ,d 2 }. Then any or- 
thogonal array with parameters OA\(n, N) over A can be used to decouple a quantum 
network consisting of n nodes of dimension d. The number of local operations used in 
this scheme is given by N . 

Proof: Let M = (m«) be an n x A r -matrix over A corresponding to the parameters 
OA(n, N). Choose unitary error bases <?i, ... , £ n where £k = {Ef : i = 1, . . . , d 2 } to 
define annihilators for each node. For each row k of M let (E^ k . : j = 1, . . . , N) be 
the corresponding local conjugations on node k. We now consider a pair of nodes, i. e., 
two rows k and I of M, and show that the local terms and the coupling between the 
two nodes are switched off. Since each pair of elements of A occurs precisely A times 
in the list ((jn,kj,mij) \ j = 1, . . . ,N) the averaged Hamiltonian H^ i on the nodes k 
and I is given by (setting Ui := E\ and Vi := E\) 

This sum is equal to zero since both annihilators are applied independently on both 
nodes. □ 
For any given number n G N of nodes there are parameters A, N such that an 
orthogonal array OA(n, N) exists. However, since we are interested in efficient schemes, 
./V has to be a polynomial in the number n of nodes. Also it is of interest to give explicit 
constructions of such schemes, i. .e, of orthogonal arrays. Whereas little is known about 
the existence of efficient schemes for general n and alphabet size s := |.A| the situation 
is much better in the case when s is a prime power. 

Corollary 1 Let an n-node quantum network with pair-interaction Hamiltonian be 
given and let the dimension d of each node be a prime power. Then there exists a de- 
coupling scheme using N local operations, where N < cn and c is a constant depending 
only on d. 



4 Note that in g the notation OA\{n, s) is used for an orthogonal array with N — As 2 in our notation, 
where s := |.4|. 
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Proof: Let s := d 2 be the size of a minimal annihilator for a d-dimensional system. 
In view of Theorem [l| we have to show that there exists an orthogonal array with 
parameters OA\(h, N) with n < h and N < cn as above. The result |], Theorem 3.20] 
gives an explicit construction for an OA((s l — l)/(s — l),s J ) for any i > 2. Hence, if 
for the number of nodes n = (s l — 1) /(s — 1) holds, we have found a decoupling scheme 
with N = s l = (s — l)n + 1 operations, i. e., N = 0(n). For general n we embed into 
an OA of this form. Switching to the next number of the form (s l — l)/(s — 1) with 
suitable i > 1 can be achieved be multiplying n with a number less or equal s, i.e., 
n < sn. □ 

Remark 1 There are tables of OAs covering the small instances (cf. ||, H; W)- We 
remark that there is a family of OAs with parameters OA(2s 1 ,2^y ~ 1) M Theorem 
6.28]. This shows that the constant c in Corollary [J can be improved to c/2. 

The following example illustrates that OAs give more efficient schemes already for 
small systems. 

Example 2 We consider the case of four three-level systems. Using the exponential 
scheme we need s 4 = 6561 local operations to decouple all interactions. Following [|| 
we obtain a decoupling scheme with the same property that uses only s 2 = 81 local 
operations. 

Selective decoupling can be achieved as follows. If the decoupling scheme is applied 
to all but one or two nodes, the remaining Hamiltonian is the local Hamiltonian of the 
node or the bipartite Hamiltonian of the two nodes, respectively. 

The assumption that every node is coupled to all the other nodes is too strong in 
many physical systems since many coupling terms might be neglected. This reduces 
the overhead for decoupling and inverting the time evolution. 

In general the interaction graph of a partially coupled network is a non-complete 
graph. Each graph can be colored, by assigning each vertex one of a number of different 
colors. Such a coloring scheme is called a proper coloring if no two connected vertices 
have the same color. The chromatic number x is the smallest number of colors required 
to properly color the graph. In a complete graph (a fully coupled network) x = n > 
but in a partially coupled network the chromatic number can be much smaller. This 
observation permits to derive more efficient decoupling schemes since if the network 
is represented by a properly colored graph, then there are no constraints on the pulse 
sequences between nodes with the same color. It is sufficient to create a decoupling 
scheme of a completely coupled x~ n ode network, and apply identical sequences to all 
nodes of the same color. 

Note that the selective decoupling scheme presented above generalize straightfor- 
wardly to the case that the dimensions of the n subsystems do not agree. Then one 
has to use so-called mixed orthogonal arrays, i.e., one has different alphabets A for 
different nodes. Although little is known about constructions of efficient mixed orthog- 
onal arrays, it is known that exponential ones exist (cf. || Section 9.3]). Furthermore, 
the mutual simulation of Hamiltonians on bipartite systems is also possible for differ- 
ent dimensions, since the proof in |2fi] does only rely on the fact that on each of the 



subsystems (possibly different) transformer groups are available. 
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2.2 Lower bound on the simulation time overhead 



In order to derive lower bounds on the overhead for the simulation time we neglect the 
free evolution and consider the weaker problem to simulate the desired Hamiltonian on 
n nodes up to local terms of each node. Note that the local terms become irrelevant 
when allowing arbitrary unitary operations on the nodes. 

In the following it will be convenient to represent the interaction Hamiltonian by 
the so-called J-matrix 



J 






Jl2 


Jl3 ■ 


• Ji n \ 


■hi 







J2n 


hi 


J32 





hn 



e M mnxmn , (6) 



\ Jnl Jn2 Jn3 / 



where the matrix J^i describes the coupling between the nodes k and I and Jik is the 
transpose of J^i in eq. ([!]). 

Let H and H be arbitrary pair-interaction Hamiltonians. We investigate the ques- 
tion whether H can be simulated by H with overhead r. This defines a quasi-order of 
the pair- interaction Hamiltonians for r = 1. A partial characterization of the quasi- 
order is expressed in terms of majorization of the spectra of the corresponding matrices 
J and J. Note that this criterion does not coincide with the criterion given in Q for 
two qubits since the latter refers to the spectrum of the Hamiltonians and not of 
the J-matrices. To make this more precise we introduce the following notation. Let 
x = (xi, . . . , Xd) and y = (yi, . . . , yd) be two d-dimensional real vectors. We introduce 
the notation j to denote the components of a vector rearranged into non-increasing 
order, so = . . . , x\), where (x[ > x<j, > . . . > x,). We say that x is majorized 
by y and write x -< y, if 

fc k 

j=l 3=1 

for k = 1, . . . , d — 1, and Y^,j=i x j = Y^j=i Vj ( see 01)- The usefulness of majorization 



techniques in quantum information processing has first been demonstrated in 1£] in 
the context of transforming states under local operations and classical communication. 
In |27]] majorization was used in the context of simulation of Hamiltonians in quantum 
networks of qubits. 

Let Spec(X) denote the spectrum of the Hermitian matrix X, i. e. the vector of 
eigenvalues. Ky Fan's maximum principle [|| gives rise to a useful constraint on the 
eigenvalues of a sum C := A + B of two Hermitian matrices: 

Spec(^ + B) -< Spec(A) + Spec(£) . (7) 

Based on these tools we derive the following bound, which generalizes the bounds given 



in [27] from qubits to qudits. The second statement of the theorem even tightens the 



bound for qubits given there. 

Theorem 3 (Lower bound) Let H and H be arbitrary pair-interaction Hamiltoni- 
ans. A necessary condition that H can be simulated with overhead r by H is that 
the spectrum of J is majorized by the spectrum of tJ. Furthermore it is necessary 
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that this majorization criterion is still satisfied after rescaling the couplings as follows: 
J' kl := SkiJki an d J' '■= SkiJ , where S = {s^i) is an arbitrary real symmetric n x n- 
matrix. 

Proof: We choose an orthonormal basis B (with respect to the trace inner product) 
for su{d). Each a £ su(d) can be represented by \a) = (tr{ao~i), . . . tr(aa d 2_i)) T with 
Oi £ B. Conjugation of a £ su(d) by a unitary u £ SU(d) corresponds to an element of 
SO{d? - 1) rotation of \a). The sub group of SO(d 2 — 1) defined by the adjoint action 
of SU(d) in this way will be denoted by 1Z. 

Representing the Hamiltonians H and H by their J matrices we see that H can 
be simulated with overhead r if and only if there is a sequence of special orthogonal 
matrices Uj = Uj\ © Ujz © • • • © Uj n , where Uj t k £ Tl for all j, k and r, > with 
J2j T j = T such that 

j = Y, t ^ ju J- ( 8 ) 



The proof now follows from Uhlmanns's theorem [21] which implies that spectrum J 
is majorized by the spectrum of r J. 

The second statement is a direct consequence from the first one since the same 
simulation procedure can be used for the rescaled problem. □ 

Note that for d > 2 not every rotation in SO(d 2 — 1) corresponds to a conjugation 
by a unitary in SU(d). Therefore this condition seems to be the weaker the greater 
d gets. A necessary condition for the fact that H can simulate H with overhead r 
is that Spec(-ff) -< rSpec(-ff), i.e. a condition based on majorization of the spectra 
of the Hamiltonians. Clearly the size of the matrices representing the Hamiltonians 
grows exponentially with n, whereas the size of their J-matrices grows only linearly 
with n. This makes the J-representation useful for calculations. Furthermore, in some 
important cases it permits to derive tight bounds (see Section || and fi~2|| ). 

2.3 Inverting 

We now consider the problem to invert an arbitrary, possibly unknown Hamiltonian in 
a quantum network, i. e. to simulate —H when H is present. It is well-known that this 
question is closely related to the construction of decoupling schemes (e. g. [22, 2(3, 17]). 



Lower bounds on the time overhead for time-reversal in n-qubit systems were given in 



1 12]. In the case of a single node we can invert the time evolution by summing over 
all elements of the unitary error basis but the identity. This trick can be generalized 
to the case of multiple nodes. For that we introduce the notion of normal form for 
decoupling schemes. 

Definition 3 (Normal form) Let O be an OA\(n,N). Assume in addition that the 
alphabet A consists of the elements of a finite group G. We say that O is in normal 
form if each entry in the first column of O is the identity element of G. 



Lemma 1 Let O be an OA\(n, N). Then there is an orthogonal array with the same 
parameters which is in normal form. 
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Proof: We identify the underlying alphabet A with an arbitrary finite group G of 
order |,4.|. Consider two rows (gi,... ,5>at) and (hi,... ,/ij\r) of O. Multiplying the 
elements of the rows by g^ 1 and h^ 1 respectively preserves the property that all pairs 
occur with frequency A since G x G is invariant under multiplication by fixed elements. 
□ 

Based on the normal form of OAs we now give an inversion scheme for a general, 
possibly unknown pair- interaction Hamiltionian. 

Theorem 4 (Inverting) Any OA(n, N) can be used to invert the time evolution of a 
quantum network consisting of n nodes. The number of local operations used is N — 1 
and the time overhead is N — 1. 

Proof: By Lemma [l] we assume that the orthogonal array O is in normal form and 
M is the corresponding n x N matrix over the alphabet A which in turn is identified 
with the elements of a finite group G. Furthermore, we choose unitary error bases 
£i, . . . ,£ n where = {Ef : % = 1, . . . , d 2 } and identify the identity element of G with 
the corresponding identity matrices. For any pair k, I of nodes we obtain the following 
identity for the averaged interaction H^ i between k and I: 

N 

~H k ,l = y^X E mk4 ® ^Vny) 1 H k,l(Em kj &> E l mij ). 

i=2 

Since the interactions between all nodes are inverted it follows that H is inverted. □ 
Similar to the results given in we derive a lower bound for the time overhead 
for inverting a dynamical evolution as follows. 

Lemma 2 (Lower bound on inverting) Let r be the greatest eigenvalue and q the 
smallest eigenvalue of the J -matrix representing H . Then fx > is a lower bound on 
the overhead for simulating —H by H . 

Proof: Denote the smallest eigenvalue of a matrix A by A m ; n (yl). Then we have 

-r = A min (-J) = Amin ( ^ Tj Uj Jiff) > rA min (J) = rq. 

j 

The inequality is due to X m i n (A+ B) > A m i n (A) + X m i n (B) for the sum of two Hermitian 
matrices A and B (see ||, Theorem III. 2). Since q is negative (J is traceless) we have 
r > —r/q. □ 
Consider a fully coupled quantum network, where each pair-interaction is of the 
form o~ a o~ a . Then the greatest and smallest eigenvalues of J are n — 1 and — 1. 
Lemma |2| gives n — 1 as lower bound on the time overhead for inversion. 

3 Harmonic oscillators with bilinear couplings 

We consider a quantum network, where the n nodes are harmonic oscillators, i. e., their 
energy values are given by E = 0, 1, 2, ... . Here we restrict our attention to the case 
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of energy values less than d and obtain as an approximation a d-dimensional Hilbert 
space for each oscillator. The Hamiltonian of the uncoupled system is given by 

n 

H :=J2hk, (9) 

k=l 

where each h k is the diagonal matrix diag(0, 1, . . . , d — 1) on the node k. Assume in 
addition that there is a coupling Hq between the individual harmonic oscillators having 
the form 

He ■= *Y^c k ia k a\ (10) 

k,l 

Here C = (cj~i) is a real symmetric n x n-matrix with zeros on the diagonal determining 
the couplings and a k and a, are the annihilation and creation operators on the k til and 
/ th oscillator, respectively. The annihilation operator a is defined by 

a|0)=0, a\E) = VE\E - 1) , E = 1, 2, . . . , d - 1 , (11) 

where {\E) | E = 0, 1, . . . , d — 1} is an eigenvector basis of the free Hamiltonian of an 
harmonic oscillator. 

Interactions of the form (|l(]) often appear if higher order terms in creation and 
annihilation operators are neglected and only that part of the total interaction term is 
considered which commutes with the uncoupled evolution corresponding to Hq. 

3.1 Decoupling with difference schemes 

The decoupling schemes presented in Section [1| can be used to decouple general pair- 
interaction Hamiltonians. In this section we consider a specific interaction which allows 
a more efficient decoupling method using the combinatorial concept of difference ma- 
trices (cf. H Section VIII]). 

The coupling He between the harmonic oscillators can be removed by transfor- 
mations of the form exp(ih k t) with t 6 R. If the time evolution according to He is 
conjugated by transformations exp(ih k t) and exp(i/i/s) on oscillator k and I, respec- 
tively, one part of the coupling term between k and / is multiplied with the factor 
exp(z(s — t)) and the adjoint part with exp(i(t — s)) since 

exjp(ihkt) exp(i/i;s) a^aj (exp(— ih^t) exp(-ihis)) = exp(i(s — t)) a k a\ , (12) 

and a similar expression holds for the adjoint term. 

We now characterize n x N matrices M having complex numbers of modulus one 
as entries which are suitable for decoupling Hq- Such a matrix M defines a sequence 
of operations on the n oscillators as follows. If (e ltl , e 4 ' 2 , . . . , e ttn ) denotes the j th 
column, this means that during the j th time interval the natural time evolution of the 
n oscillators is conjugated by the local transformation 

n 

Y[ exp(ih k t k ) . (13) 
fc=i 
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The total effect of this scheme is that the term akaj obtains the factor (rrty|m&), 
where is the k th row of M and (-|-) is the usual inner product in C^. This gives 
the following decoupling criterion: all couplings are removed if and only if the rows 
are orthogonal in the usual sense. There is a canonical way of finding n vectors having 
this property by taking the Fourier transform of the standard basis of C n . However, 
the rotations exp(2irihj /n) required to be implemented are very close to the identity 
for large n. 

An alternative way of constructing such a matrix is given by difference schemes 
ID, |9|. We choose the numbers tj, in each row to be of the form 2irr/u where u > 2 is a 
natural number and r is an element of the cyclic group Z u = {0,1,... , u — 1} of order 
u. 

Definition 4 An nx N array D(n, N) with entries in Z u is called a difference scheme 
based on Z u if the difference vector of any two rows has the property that each element 
of Z u occurs equally often. 

Let D(n,N) be a difference scheme based on Z u . We construct annxJV complex 
matrix M from D by replacing each entry r by e 2mr ^ u . The rows of M are vectors in 
C N and since D is a difference scheme they are orthogonal, i. e., 

MM* = N1 N . (14) 



Therefore, M satisfies the decoupling criterion. In view of (14), a difference scheme 
D(n, N) in which n = N is also called a generalized Hadamard matrix |^, Section VIII] 
of order n over Z u . In particularly, any (ordinary) Hadamard matrix of order n is a 
difference scheme D(n,n) over Z2. 



3.2 Recoupling disjoint cliques 

The scheme presented above does not only allow to remove all interactions but also to 
achieve the following selective recoupling without time overhead. Partition the set of 
nodes into n' disjoint subsets (called "cliques") and remove only the couplings between 
nodes in different cliques. This can be achieved by applying the same sequences of 
transformations on all nodes in the same clique, since this does not affect the interac- 
tions among them. Then it is sufficient to construct a difference scheme with only n' 
rows since each row refers to one of the cliques. Note that an analogous way of "clique 
decoupling" is also possible for the following kind of n-qubit interaction. Assume that 
all qubits are coupled by the interaction a x (g> a x + a y (g> a y + a z <g> o z . Then the interact- 
ing is invariant with respect to simultaneous unitary rotations on both qubits. Hence 
decoupling schemes for n' qubits define a "clique decoupling" scheme for n' cliques. 



3.3 Simulation of different coupling strengths 

An interesting problem arises if we want to switch off interactions between arbitrary 
pairs or, even more general, to weaken the interaction between some of the oscilla- 
tors. These goals are special instances of the general problem to simulate a coupling 
Hamiltonian using Hq where C is an arbitrary real symmetric coupling matrix. 
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Now we define a matrix T in such a way that the entry-wise product (Schur-product) 
of T and C is C. Of course this is not possible if an entry of C is equal to zero and 
the corresponding entry of C does not vanish. This corresponds to the fact that one 
cannot simulate a coupling between nodes which are not coupled. We choose the vectors 



irik introduced in Subsection 3.1 with complex numbers of modulus one such that the 
corresponding Gram matrix G = {{mk\mi))k,i=i,... n coincides with T on all off-diagonal 
entries. Note that the diagonal entries of G always give the number of time steps of 
the simulation. One can generalize this by choosing time steps of different length and 
define the vectors m k as mappings from the interval [0, t] to the complex numbers of 
modulus one. By taking the inner product 

(m k \mi) := / m k {t')mi(t')dt' 



JO 

we obtain the same statement as above: the time overhead t is determined by the 
diagonal entries of G. Note that in our formulation a weak interaction can be used to 
simulate a strong one since we allow time overhead for the simulation. 

A lower bound for this overhead is given by the absolute value of the least eigenvalue 
of T since G being a Gram matrix is positive. This coincides with the eigenvalue 
criterions stated in [27| derived for n spin 1/2 systems. The lower bound for simulating 



its own inverse has been shown to be n — 1 for the spin-spin interaction of the form 
o~ z <S> o~ z [27, [l^j . In strong analogy, we obtain the same bound n — 1 for simulating 



—Hq = H-c by Hq when C has vanishing terms only on the diagonal, since in 
this case the corresponding matrix T has only entries —1 except for the diagonal 
entries which are zero. This proves the bound since the absolute value of the least 
eigenvalue of T is n — 1. Note that lower bounds based on the T-matrix do only refer 
to simulations using unitary transformations of the form exp(ihjt) with t £ R. The 
following theorem shows that the lower bound n — 1 for time reversal is even valid if 
general unitary transformations on each oscillator and presents an inversion scheme 
attaining this bound. 

Theorem 5 (Optimal inversion) A time optimal scheme for inverting a coupling 
of the form He as in equation is given by the vectors 



mk ._ t^lmk/n ^ e 2m2k/n^ ^ e 2m(n-l)/n\ 

Proof: First note that due to the length of the vectors this simulation has time 
overhead n — 1. In order to prove that this is optimal, define the following linearly 
independent elements of su(d). 

X r := \r)(r - 1| + \r - l)(r| 

and 

Y r := i\r)(r — 1| — i\r — l)(r| 

for r = 1, . . . ,d— 1. These matrices are orthonormal with respect to the inner product 
(V|W) := tr(VW)/2. One may supplement these 2d — 2 vectors to a orthonormal 
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basis of su(d), but the completion is irrelevant since the interaction among each pair 
of oscillators can be written as an expression in X r and Y r . 



a (g> a+ + a) ® a = V^(^r ® X s + T r (g> F s ) . 

The coupling matrix J (see eq. (@)) can be constructed as follows. Define a {2d — 2) x 
(2d - 2)-matrix A by 

A := |0>^| + |^)^|, 

with 

|0) := (vT, V2,...,Vd- 1,0,...,0) T 

and 

|V) := (0,...,0,VT,^,...,v / d^T) T . 

With respect to the basis described above, for all pairs (k, I) of oscillators the coupling 
matrices J^i are the same and given by completing A by embedding into a (d 2 — 1) x 
{d 2 — l)-matrix A'. Since the vectors \4>) and are orthogonal, the spectrum of Jk,i 
is a, a, 0, ... ,0, where a = {(fr\4>) = > 0. The spectrum of J contains two copies 

of the value (n — l)a and n — 1 times the value —a. The other eigenvalues are zero. 
This can be seen by writing J as a tensor product M <g> A' where M is an n x n-matrix 
with entries on the diagonal and 1 elsewhere. The smallest value r > such that 
the spectrum of r J majorizes — J is n — 1. This proves optimality of the time reversal 
scheme. □ 
An upper bound for the general simulation problem can also be derived in strong 



analogy to 27] and is given by the so-called weighted chromatic index Wt of the matrix 



T. This concept has been introduced in [11] in a related context in order to quantify 



the complexity of a general pair-interaction Hamiltonian on n qubits. For each s > 
define a graph G s on the n nodes as vertices which has an edge (k, I) if and only if the 
absolute value of the entry T^i is greater than s. Let n s be the chromatic number of 
G s (see pi), i.e., the number of colors required for coloring the edges of G s in such a 
way that no two edges with a common node receive the same color. Then define Wt 
as 



oo 



Wt ■= / n s ds 



o 

which defines a generalization of the chromatic index for weighted graphs. The key 
idea to prove this upper bound can easily be understood if one assumes T to have 
entries of modulus one or zero. Then Wt is the chromatic index of a graph indicating 
which couplings should not be removed. Given an admissible coloring of this graph, 
we define a simulation with Wt steps as follows. Each color c defines a step in which 
we remove all those interactions which are not colored by c. This step can be executed 
without time overhead as explained in Subsection [T^. Hence in each step, only the 
couplings between disjoint oscillator pairs remain. 
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4 Comparison with other methods 



In this section we relate the decoupling method based on orthogonal arrays which 
was presented in Section § to the approach of [17| for the case of decoupling in qubit 



networks (see also [24] for decoupling in qubit networks based on orthogonal arrays). 
Accordingly, let the system Hamiltonian be written in the form 



kl;a/3 



(15) 



k:a 



where a a are the Pauli matrices. 

We explain briefly the approach of [17] to construct decoupling schemes. In each 
interval, each a\ acquires a + or — sign, which is controlled by the applied local 
unitaries to be described. The coupling Jki-a^a^p for k is unchanged (negated) if 
the signs of cr^ and cr^ agree (disagree) . Note that the signs of the three Pauli matrices 

acting on the same qubit k are not independent. 

Conjugating with the transformations l k ,a x ,a^, a l x the acquired signs for a x , o~y,a z 
are given by (+ + +), (+ - -), (- + -), (- - +). 

Following 0] a decoupling scheme for n qubits that concatenates N intervals can 
be specified by three nx N sign matrices S x , S y ,S z , related by the entry- wise product 
S x * S y = S z . We say that the three matrices satisfy the Schur condition since this 
entry-wise product is usually called the Schur product. The (k,j) entry of S a is the 
sign of o~ k in the j th time interval. Therefore decoupling is achieved if any two rows 
taken from S x , S y , S z are orthogonal. 

The following theorem establishes a connection between decoupling schemes con- 
structed using orthogonal arrays as described in Section |2| and the decoupling schemes 
specified by sign matrices S x , Sy, S z . 

Theorem 6 A decoupling scheme constructed using an orthogonal array OA(n, N) 
over the alphabet A = {1,2,3,4} and the Pauli basis £ = {1, cr x , o~ y , a z } for all nodes 
gives rise to sign matrices S x ,S y ,S z satisfying the Schur and orthogonality conditions. 



Proof: We identify the operators of £ with the elements of A according to 1 i— ► 1, 
2 i— > <r x , 3 i— ► a y , and 4 i— ► a z . Conjugating with the operators of £ the Pauli matrices 
acquire the following signs: 





1 


2 


CO 


4 


0~ x 


+ 


+ 






Oy 


+ 




+ 






+ 






+ 



(16) 



Starting from the given orthogonal array we now construct the three sign matrices 
S x , Sy, and S z . Pick any two rows k and I of the OA. We may assume that the two 
rows have the following form (or else we apply a suitable permutation of the columns) 



1 1 


1 


1 


2 


2 


2 


2 


CO 


1 2 


CO 


4 


1 


2 


3 


4 


1 



















3 
2 

A times 



4 4 4 4 
12 3 4 



(17) 
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since all pairs appear equally often (A times) in the OA. Let A = (+ + •••+) be the 
vector of length A = N/16 containing only +. By substituting the entries of the rows 



k and the I of the OA by the corresponding sign assignments in Table (16) we define 
the following six rows of S x ,S y ,S z , respectively: 



Sx;k 


= A (8) 


(+ + — ) 


<8 


(+ + ++) 




= A§ 


(+ - +-) 


<8 


(+ + ++) 


Sz;k 


= A® 


(+ - -+) 


<8 


(+ + ++) 


S x ;l 


= A <8 


(+ + ++) 


® 


(+ + — ) 


Sy-l 


= A (8 


(+ + ++) 


<8 


(+"+-) 


S Z ;l 


= A® 


(+ + ++) 


<8 


(+ - "+) 



Obviously, S x , S y , and S z are orthogonal and satisfy the Schur condition S x * S y = S z . 
□ 

Finally, we give an alternative proof for the existence of a decoupling scheme for 
n = (2 2m — l)/3 qubits using N = 2 2m time intervals. A decoupling scheme with these 
parameters can be constructed using orthogonal arrays [^4| and Hadamard matrices 

Let V be the vector space F™, where rn > 1 and let F4 = {0, 1, u, to 2 = l+u}, where 
= 1, be the Galois field with 4 elements. Recall that the number of d-dimensional 
subspaces of an m-dimensional vector space over ¥ q is given by 

{q m - l)^™- 1 -!)••• ( q m ~ d + l - 1) 



m 
d 



(q d -l){q° 



1) 



(9-1) 



(18) 



(cf. & Lemma 2.14, Section I]). For the special case q = 4 and d = 1 formula (fill) 
shows that there are (4 m - l)/(4 - 1) = (2 2m - l)/3 lines in Wf. Note that different 
lines intersect in the point {0} only. Hence by taking the set of all one-dimensional 
subspaces of F™ we obtain a maximal spread in F™, i.e. a collection of subspaces Ui 
partitioning F™ with the additional property that 

Ui n Uj = {0} . 

We define a map ip from F4 onto { — 1, +1} 4 as follows: 



p(0) 

This is the Hadamard matrix H4 = 



(+1,+1,+1,+1) 
(+1,-1, +1,-1) 
(+1, +1,-1,-1) 
(+1,-1,-1, +1) 



Ho 



H2 where H2 is the usual Hadamard matrix 
of size 2. Therefore all rows are orthogonal]^ Note that the last three rows satisfy the 
Schur condition. 



The fact that this matrix is indeed the Hadamard matrix can also be derived with the help of group 
characters. More precisely, we consider F4 as a two-dimensional vector space over F2 and let tr denote the 
trace map of this field extension |l0|, Section 4.15]. For all z 6 F4 the map <p(z) : x 1— > (— l) tr ( za: ) is an 
irreducible character of the additive group (F4,-)-) (which is isomorphic to Z2 x Z2). Hence, orthogonality 
of the rows follows from the orthogonality of the characters. 



16 



We extend the map (p to vectors v = . . . , v m ) S FJ 1 by denning the map 

<j>(v) := pfa) ® • • • ® <p(v m ) G {-1, l} 4m . (19) 

The image of is the set of all rows of the Hadamard matrix Hf 2m . Let U k = (v k ) be 
a maximal spread of F™. By evaluating <p on the three elements of U k (except for the 
zero vector 0) we get three orthogonal vectors satisfying the Schur condition. We can 
take them as rows of S x , S y , S z : 



S X ;k — 


■ V k ) 


Sy\k = 


<p(u; 2 • v k 


S Z ;k = 


0(1 • v k ) 



This shows that the second through the last rows of the Hadamard matrix H® m can 
be divided into (2 2m — l)/3 disjoint 3-subsets, each with rows that satisfy the Schur 
condition. The rows in a 3-subset can be chosen as rows of S x ,S y and S z , respectively. 

5 Conclusions 

We have shown that pair-interactions between the subsystems of a multipartite quan- 
tum system can be decoupled efficiently if a sufficiently large set of local control op- 
erations on the subsystems is available. Such decoupling schemes can be constructed 
using orthogonal arrays (a concept of combinatorics). The rows of these arrays de- 
fine pulse sequences of local operations taken from a unitary error basis. We discuss 
the connection between the decoupling method based on orthogonal arrays and those 
introduced in [|l7]]. 

We have shown that mutual simulation of pair-interaction Hamiltonians in multi 
node systems is possible provided that a so-called transformer group of transformations 
is available.^ The upper bound 0{n 2 ) on the simulation time is a consequence of 
the existence of selective decoupling schemes. The construction of a time-optimal 
simulation leads to a non-trivial convex optimization problem. We have derived a 
lower bound on the time overhead in terms of the spectrum of the matrix describing the 
coupling parameters. For some interactions simpler decoupling schemes can be devised: 
for bilinear coupling of harmonic oscillators selective decoupling can be achieved using 
so-called difference schemes. The condition for time optimality of mutual simulation of 
different bilinear couplings can concisely be expressed in terms of linear algebra. Based 
on this we have constructed time optimal schemes for time reversal. 

From the results shown in this paper it follows that the time optimal implementa- 
tion of unitary transformations turns out not to be a matter of optimal factorization 
into parallelized bilocal quantum gates alone. The transformation has rather to be 
written as the solution of a time-dependent Schrddinger equation where the occurring 
Hamiltonians are those which can be simulated with small time overhead. This leads 
to another definition of quantum complexity different from the discrete one measured 
by counting the number of elementary gates. 

6 Note that similar results have been developed in p0| , independently. 
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